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TEXTURE FUNCTIONS IN IMAGE ANALYSIS; 

A COMPUTATIONALLY EFFICIENT SOLUTION 

ABSTRACT 

A set of statistics that measures^visually perceivable! textures in images by use of co-occurrence 
matrices has previously been developed. Presented here is a computationally efficient means for 
calculating texture measurements from digital images by use of the co-occurrence technique. This 
paper discusses the calculation of the statistical descriptors of image texture and presents a solution 
that circumvents the need for calculating and storing a co-occurrence matrix. The results show that 
existing efficient algorithms for calculating sums, sums of squares, and cross products can be used to 
compute complex co-occurrence relationships directly from the digital image input. 
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TEXTURE FUNCTIONS IN IMAGE ANALYSIS: 

A COMPUTATIONALLY EFFICIENT SOLUTION 

1.0 INTRODUCTION 

This paper presents a computationally efficient solution developed as part of the Multispectral 
Linear Array Supporting Science Studies (MLASSS) at the Goddard Space Flight Center (GSFC), 
for calculating texture statistics from digital images. Presented here are the mathematical founda- 
tions of image-texture calculations based on the Spatial Grey Tone Dependence (SGTD) method 
developed by Haralick (1) and Haralick et al. (2) and reviewed by Conners (3) and Conners and 
Harlow (4). 

One aim of the MLASSS is to determine the synergistic effects of increased sensor spatial, 
spectral, and radiometric resolution. Spatial resolution studies at GSFC have focused in increases 
in image information content with increased spatial resolution and evaluation of sensor systems 
with mixed spatial resolution as possible candidates for future land remote-sensing missions. To 
this end texture analysis, in particular SGTD, was seen as one way to quantify the increased spatial 
information apparent in high-resolution digital imagery. This paper documents the development 
and use of software at GSFC to implement Haralick’s algorithms. 

2.0 BACKGROUND 

Landsat-4’s successful launch and subsequent successful operation of the Thematic Mapper 
(TM) and Multispectral Scanner (MSS) heralded a new era in space-based remote sensing. Advances 
include improved spatial resolution of 30-meters instantaneous field-of-view, 7 spectral bands and 
8-bit quantization for the TM. In the future, multispectral linear array (MLA) technology will 
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make it possible to obtain digital imagery of vastly improved spatial resolution, 10 to 15-m in the 
visible and near infrared (5). In anticipation of this increased capability, several investigators have 
used high-resolution aircraft scanner data to study the tradeoffs associated with increased spatial 
resolution, processing strategies, and costs. A potentially serious problem has been encountered in 
the use of conventional unsupervised or supervised per-pixel classifiers: as spatial resolution in- 
creases, classification accuracies tend to diminish in areas of high spatial complexity (6). Further- 
more, as the proportion of mixed pixels increases, classification accuracies decrease. An example 
would be the decreasing classification accuracy of small agricultural fields whose size approaches 
the sensor IFOV. Conversely, heterogeneous land covers characterized by small high-frequency 
components tend to be averaged at lower resolutions so that classification accuracies are higher with 
per-pixel classifiers. These results have been independently confirmed by Latty (7) for forested 
sites. 

These studies point to a need to incorporate spatial information in the classification process. 
Several methods have been advocated in addition to Haralick’s SGTD method; they include spatial/ 
spectral context, used by Tilton and Swain (8), and categorical/spatial context, developed by 
Wharton (9). This paper discusses the Haralick SGTD algorithm and derives a computationally 
efficient means for calculating various texture statistics derived from a spatial gray-tone co- 
occurrence matrix. 

Texture analysis as discussed in this paper is used to quantify the spatial information in a 
digital image by measuring the spatial arrangement of gray tones within it. The recent literature 
includes a review of various texture-analysis methods by Haralick (2) and an update by Davis (10). 
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Conners and Harlow ,(4) investigated the theoretical merits of various texture-analysis strategies 
for quantifying image patterning. Cox et al. (1 1) and Weszka et al, (12) conducted empirical com- 
parisons of various texture measures. 

The SGTD method has been used frequently by investigators working with remotely sensed 
data, including Haralick et al. (1), Hsu (13), Jensen (14) and Toll (15), Schowengerdt (16), and 
Weszka et al. (12). Compared to first-order statistics such as mean and standard deviation, SGTD 
has greater potential but is computationally more complex. The SGTD method transforms the 
gray values within a neighborhood (window) into a two-dimensional gray-tone co-occurrence of 
gray-tone pairs i and j as measured among angle a for distance d and can be interpreted as a prob- 
ability matrix of gray-tone pairs. Haralick et al. (1) introduced a number of statistics based on 
information theory to describe such matrices, and as part of the spatial studies of the MLASSS, 
eight have been put to use at GSFC on the HP-3000-based Interactive Digital Image Manipulation 
System (IDIMS) in the Applications Directorate. What follows is a discussion of these eight algo- 
rithms and their implementation. 

3.0 CO-OCCURRENCE CALCULATIONS 

Given a rectangular matrix (window) of values (brightness), an occurrence matrix is defined 
as the frequency with which a value of i precedes a value of j in the direction a. Call this matrix 
qij(«). 

A co-occurrence matrix is defined as 

Pjj(a) = qjj(a) + qy(cdvr) 

or the sum of the occurrence matrix in one direction plus the occurrence matrix in the opposite 


direction. 
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It should be noted that 
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qy(orHr) * ^(a) 

and * 

qjj(a) = qjiCa-hr) 

That is, the number of times i precedes j in direction a is exactly the number of times j pre- 
cedes i in the opposite (ortir) direction. 

Consequently 

= qy(a) + qji(a) 

or alternatively 

p „(a) = qy (a) + q'ij(«) 

Where q*j(a) is the transpose of q^(a). 

Further relationships between the co-occurrence matrix and the occurrence, or precedence, 
matrix can be readily seen. 

The sum of all elements within the matrix 

NI NJ NI NJ N1 NJ 

Z S p ij(a) = Z Z qjj(a)+ Z Zqjj(a) 

i j ‘ j i i 

NI NJ 

-2 Z Z qjjCa); due to symmetry 
‘ j V 

Furthermore, since the precedence matrix (q) is strictly a count of relationships within the 
window of precedence, the sum of that matrix is simply the number of relationships within that 
window. Since each cell specifies a single relationship within the precedence matrix 


4 




original page is 

OF POOR QUALITY 


Nl NJ 

I ZPnCa) * 2N a ; where N~ is the number of a-pairs 
l j u 


In a similar fashion 


and 


sa 1 *p # -iz,*q # + sa 1 *qj, 

S2 u p u 55 + S2 ji^U 

- 2 SX^q^; due to symmetry. 


The occurrence matrix q is more literally a precedence matrix in that it is a count of the num» 
ber of times a value i precedes a value j in direction a. The particular values that precede occur in 
a subset of the original window, depending on the direction a. Specifically, in the diagrams below 
the range of precedence is indicated. 



0° 90° 45° 135° 

Range of Precedence 

Similarly, the range of succession depends on the direction a. 



0 ° 90° 45° 135° 


Range of Succession 
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It should be noted in these diagrams that the number of rows (scan lines) and the number of 
columns (pixels) within the range of precedence both depend on the direction a. However, for any 
given direction a the limits on the range of precedent are equivalent to the limits on the range 
of succession. 

The average brightness value within the range of precedence can be expressed in two* ways. 
The brightness value times the relative frequency of occurrence of that brightness-summed-is one 
way to express the mean: 

nv 

M ■ Z tfi 

i«l 

Expressed in the language of precedence and co-occurrence, that relative frequency is the 
relative number of times which that brightness (i) preceded all other brightness values. 

NV 

Z qij 

fi = — where N~ is the total number of a=pairs 

iNj r* ** 


Thus, one expression for the mean is 

A more efficient method of calculating the average brightness within the range of precedence 
is only to sum all brightness values and divide by the number of a-pairs. 


NS NR 


M Z Z x kl/^a 

k 1 


Where NS and NR are the limits of the range of precedence. 
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The equivalence of these two measures of average brightness points ouS that 

NV NV NS NR 

1 21 irljj = ^ ^ ^kl 

l j w k > 

In a similar fashion it can be demonstrated that for the range of succession 

NV NV MS MR 

ZZM^ Z Z X k , 

* I j * 1 

Where MS and MR are the limits of the range of succession. 

Again, for higher order moments of precedence or succession it can be shown that 

S2i^ = 22Xj‘ 1 

and 

SStiqy = Sr k iX k+ail+b 
where a and b depend on the direction a. 

Sum, sums of squares, and crossproducts are computationally efficient algorithms. These 
translations have made it possible to compute fairly complex co-occurrence relationships by use of 
computationally efficient techniques. 

4.0 CALCULATION OF TEXTURE STATISTICS 

Eight texture functions based on the pairs of co-occurrence were developed. The background 

and application of these functions has been fully described elsewhere (1 , 2, 3, 4). 

Of the eight functions, three (variance, skewness, kurtosis) can be phrased readily in a com- 
putationally efficient form. Three further functions (difference moment, homogeneity, correla- 
tion), which are based on the co-occurTence matrix approach, can be reduced to a more efficient 
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and standard form. The last two factions (energy, entropy) have not yielded an efficient solution 
and still require co-occurrence calculations. 


4.1 Variance, Skewness, Kurtosis 

Within each window the first four moments of the brightness values were calculated. Define 


Mr 


>-zx r 


Then 


Mi - Mi 

Ma - Mi - (Mi) 2 

Ms * Ma *“ M 1 + 2 (m1 ) 3 

M4 * Mi “ 4jLt 3 Mi + 6 m 2 (M i ) 2 ” 30i ) 4 


Using those definitions, the variance is 




(N-l) 


The coefficient of skewness is defined as 


/r — _ M3 

(M 2 ) 3/2 


And the coefficient of kurtosis is given as 


h -if* 


These statistics require only that the four raw moments (fi' r ) be accumulated for each pixel 
within a window as a whole. 
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The co-occurrence formulation of the difference moment function is 

NV NV 

2,2, (H) 2p ij 


Since this is a symmetric function [(i-j) 2 * (j-i) 2 1 occurrence, or precedence, matrix formu- 
lation is 


NV NV 

0 - 2 2 I (i-j) 2 qy 

i j y 

This is again the standard grouped-data formulation of a moment. Writing it in more efficient 
ungrouped format, 

SS(X.-X a ) 3 

N a 

Where the X t values are the brightness values of preceding pixels and the X 2 values are the 
values of their successors. The summation is over the full range of precedence. 


4.3 Homogeneity 

The co-occurrence formulation of homogeneity is 


NV NV 

r= I I 


i 


j 1 + (H) 


-H2 U 


Once again we can observe the symmetry, and the precedence formulation is similar: 

NV NV * 

T = 2 I I 1 + (i-j) 2 q« 

i j 

And finally, tiu more standard ungrouped form of this moment equation: 
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4.4 Correlation 

The co-occurrence formulation of this function has been expressed as 


^ Z (i-m) (j-m) Py 

P Z Z (i-m) 2 Py 

i j u 

where m = 22ipy/22py 
In expanded form 

ISijpy - mSSipy - mZSjpy + m 2 22py 
P £21* py - 2mSSipy + m 2 22py 

22jpy = SSipy since Py is symmetric 

and SLipy * m2 Spy 
therefore 

SSaPij-m^Spjj 

SSi’py-m'SSpy 

In terms of the precedence matrices the correlation can be expressed 

222ijpy-2m 2 N a 
' 5 ~22i J q ij + IIrq ij , -2m 2 N a; 

And in ungrouped terms 
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In addition 


ssx , x 2 


N, 


- m 2 


a 


1JB2X? 22X| x . 

.(-T~ +-TT- i )-m- 


T N, 


a 


N 


a 
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m 8 

2N a 


(■ 


SSiq^j + 22jqjj > , 


2N 


) 3 


a 


SSXj, ZSX 3 
2N„ 2N a 


which is the square of the average of the mean precedence and the mean successor value. 


4.5 Energy and Entropy 

Neither of these functions has been found amenable to any simplification. The energy func- 
tion has been expressed 


and the entropy function 


NV NV 

E t - Z Z Py' 


NV NV 

= zs 

i i 


[-Pij log Pjj 


As can be seen, both functions use nonlinear forms of the frequencies (not of the brightness 
values). This nonlinearity requires that the entire co-occurrence matrix be accumulated before the 
function can be evaluated. Since there is little limitation on the size of NV (spectral variance), 
these functions are notably expensive in terms of computer resources. 
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5.0 SUMMARY 


We have presented a straightforward method of computing texture statistics from digital 
images. Based on the SGTD method, the computation of co-occurrence matrices is viewed as 
matrices of precedence and succession from which spatial/spectral relationships between neighbor- 
ing pixels can be calculated. The precedence/succession method allows for the direct calculation of 
several texture statistics directly from the input data window without the need to calculate a time 
co-occurrence matrix as an intermediate step. This allows the realization of certain economies in 
memory and processing speed through the elimination of the co-occurrence matrix in the computa- 
tion. However, instances in which individual co-occurrence matrix elements must be operated 
upon, the precedence/succession technique is not applicable. 
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